#ifndef _RPA_H_
#define _RPA_H_
#include<complex>

#include"system_hf.h"
#include<eigen3/Eigen/Eigenvalues>

using std::max;
using std::min;

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;
using Eigen::ComplexEigenSolver;
using Eigen::MatrixXd;
using Eigen::MatrixXcd;
using Eigen::VectorXd;
using Eigen::VectorXcd;

typedef complex<double> complexd;

///< particle and hole twobody states, particle first
class TwoBodyState_ph: public TwoBodyState
{
 public:
  TwoBodyState_ph(int p,int h):TwoBodyState(p,h){}
  //below method is not used in rpa cal.
  bool isInAChannelWith(const TwoBodyState &) const{}
  virtual bool operator<(const TwoBodyState&) const{}
};


class RPA
{
 public:
  typedef vector<System_Tz_HF::OrbitalType> HFOrbitals;
  
  RPA(const System_Tz_HF* Ptr);

  int set_ph_same(int J,int par,vector<TwoBodyState_ph>&ph);
  void set_AB_same(int J,int par);//set matrix A and B which particle and hole are same kind of particles


  int set_ph_np(int J,int par,vector<TwoBodyState_ph>&ph);
  int set_ph_pn(int J,int par,vector<TwoBodyState_ph>&ph);
  void set_AB(int J,int par);//set matrix A and B which particle and hole are different kind(neutron,proton), charge exchange procedure.

  void setTDAflag();//for TDA cal., set B=0 in rpa eq.
  bool isTDA;
  
  void cal_same(int J,int par);///<par=0 for +, par=1 for -
  void cal_ce(int J,int par);///<par=0 for +, par=1 for -
  void cal(int type,int J,int par);///<type=1, call cal_ce; type=0, call cal_same
  void cal();///< cal all

  double getDeltaE();///< return correlation energy.
  double calOccuNum();///< cal. occupation number of orbitals. Phys. Rev. 175,4(1968)
  vector<double> OccuNum;

  void set_q(int J,int par,int T);///<set up onebody matrix element for Q;
  double BEJ(int J,int par,int v);///<B(EJT,0->v)
  double S(int J,int par,int T);
  void R(int J,int par,int T,ostream&,const double gamma=2,int points=600,const double dE=0.1);
  void set_q_for_rou(int J,int par,int tz,double r);
  double rou_tr(int J,int par,int tz,int v,double r);///< transition density
  double print_rou_tr(int J,int Par,double omega,int points,double dr=0.1);
  vector<double> q;///<onebody matrix element for Q;

  void set_q_for_Beta_F(int type);///< set up onebody matrix ele. for Fermi beta decay. type=1 for plus. type=-1 for minus
  void set_q_for_Beta_GT(int type);///< set up onebody matrix ele. for GT beta decay.

  //reduced matrix element
  double B_F(int type,int v);/// 0+ -> 0+
  double B_GT(int type,int v);///0+ -> 1+
  double S_F();//cal. sum rule for Fermi
  double S_GT();//cal. sum rule for GT

  
  void print_detail_info(int J,int Par,int T);
  void print_XY(int J,int Par,int v);

  void print_detail_info_ce(int J,int Par,int Tz);
  void print_XY_ce(int J,int Par,int Tz,int v);
  
  vector<vector<TwoBodyState_ph> >ph_sames;//both the particle and hole are protons or neutrons
  vector<vector<TwoBodyState_ph> >ph_nps;//neutron particle and proton hole
  vector<vector<TwoBodyState_ph> >ph_pns;//proton particle and neutron hole

  
  const System_Tz_HF* pSystem_Tz_HF;
  vector<MatrixXd> A_sames;
  vector<MatrixXd> B_sames;
  vector<MatrixXd> X_sames;
  vector<MatrixXd> Y_sames;
  vector<VectorXd> Omega_sames;


  
  vector<MatrixXd> Anpnps;
  vector<MatrixXd> Apnpns;
  vector<MatrixXd> Bnppns;
  
  vector<MatrixXd> T_plus_Xs;
  vector<MatrixXd> T_plus_Ys;
  vector<VectorXd> T_plus_Omegas;

  vector<MatrixXd> T_minus_Xs;
  vector<MatrixXd> T_minus_Ys;
  vector<VectorXd> T_minus_Omegas;

  int Jmax;  
 private:
  double a(int bra_p,int bra_h,int ket_p,int ket_h,int J);
  double b(int bra_p,int bra_h,int ket_p,int ket_h,int J);
};
#endif
